#include "euler.hh"
#include <map>
#include <math.h>
namespace Dune {
  namespace Equation {
    Euler_<1>::Euler_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    Euler_<2>::Euler_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    Euler_<3>::Euler_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    Euler_<1>::~Euler_(){
    }
    Euler_<2>::~Euler_(){
    }
    Euler_<3>::~Euler_(){
    }

    void Euler_<1>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("momentum1");
      varname.push_back("energy");  
      set_equ_para("gamma",5.0/3.0);
      para["gamma"]=ptree.get("equation.gamma",5.0/3.0);
    }
    void Euler_<2>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("momentum1");
      varname.push_back("momentum2");
      varname.push_back("energy");  
      set_equ_para("gamma",5.0/3.0);
      para["gamma"]=ptree.get("equation.gamma",5.0/3.0);
    }
    void Euler_<3>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("momentum1");
      varname.push_back("momentum2");
      varname.push_back("momentum3");
      varname.push_back("energy");  
      set_equ_para("gamma",5.0/3.0);
      para["gamma"]=ptree.get("equation.gamma",5.0/3.0);
    }

    void Euler_<1>::get_equation_flux(int whichdim, double* fequ, double* var){
      get_1d_fx(fequ,var);
    }
    void Euler_<2>::get_equation_flux(int whichdim, double* fequ, double* var){
      if(whichdim==0)
        get_2d_fx(fequ,var);
      if(whichdim==1)
        get_2d_fy(fequ,var);
    }
    void Euler_<3>::get_equation_flux(int whichdim, double* fequ, double* var){
      if(whichdim==0)
        get_3d_fx(fequ,var);
      if(whichdim==1)
        get_3d_fy(fequ,var);
      if(whichdim==2)
        get_3d_fz(fequ,var);
    }
    
    double Euler_<1>::get_thermal_pressure(double* var){
      double pressure;
      pressure=(para["gamma"]-1.0)*(var[2]
      -1.0/2.0*std::pow(var[1],2)/var[0]);
    }
    double Euler_<2>::get_thermal_pressure(double* var){
      double pressure;
      pressure=(para["gamma"]-1.0)*(var[3]
      -1.0/2.0*std::pow(var[1],2)/var[0]
      -1.0/2.0*std::pow(var[2],2)/var[0]);
      return pressure;
    }
    double Euler_<3>::get_thermal_pressure(double* var){
      double pressure;
      pressure=(para["gamma"]-1.0)*(var[4]
      -1.0/2.0*std::pow(var[1],2)/var[0]
      -1.0/2.0*std::pow(var[2],2)/var[0]
      -1.0/2.0*std::pow(var[3],2)/var[0]);
      return pressure;
    }

    double Euler_<1>::local_max_eigenvalue(double* var){
      double pressure = get_thermal_pressure(var);
      return std::sqrt(para["gamma"]*pressure/var[0]);  
    }
    double Euler_<2>::local_max_eigenvalue(double* var){
      double pressure = get_thermal_pressure(var);
      return std::sqrt(para["gamma"]*pressure/var[0]);  
    }
    double Euler_<3>::local_max_eigenvalue(double* var){
      double pressure = get_thermal_pressure(var);
      return std::sqrt(para["gamma"]*pressure/var[0]);  
    }
    
    void Euler_<1>::get_1d_fx(double* fequ, double* var){
      double pressure;
      pressure=get_thermal_pressure(var);
      fequ[0]=var[1];
      fequ[1]=std::pow(var[1],2)/var[0]+pressure;
      fequ[2]=(var[2]+pressure)*var[1]/var[0];
    }
    
    void Euler_<2>::get_2d_fx(double* fequ, double* var){
      double pressure;
      pressure=get_thermal_pressure(var);  
      fequ[0]=var[1];
      fequ[1]=std::pow(var[1],2)/var[0]+pressure;
      fequ[2]=var[1]*var[2]/var[0];
      fequ[3]=(var[3]+pressure)*var[1]/var[0];      
    }
    
    void Euler_<2>::get_2d_fy(double* fequ, double* var){
      double pressure;
      pressure=get_thermal_pressure(var);    
      fequ[0]=var[2];
      fequ[1]=var[2]*var[1]/var[0];
      fequ[2]=std::pow(var[2],2)/var[0]+pressure;
      fequ[3]=(var[3]+pressure)*var[2]/var[0];  
    }
    
    void Euler_<3>::get_3d_fx(double* fequ, double* var){
      double pressure;
      pressure=get_thermal_pressure(var);    
      fequ[0]=var[1];
      fequ[1]=std::pow(var[1],2)/var[0]+pressure;
      fequ[2]=var[1]*var[2]/var[0];
      fequ[3]=var[1]*var[3]/var[0];    
      fequ[4]=(var[4]+pressure)*var[1]/var[0];
    }
    
    void Euler_<3>::get_3d_fy(double* fequ, double* var){
      double pressure;
      pressure=get_thermal_pressure(var);      
      fequ[0]=var[2];
      fequ[1]=var[2]*var[1]/var[0];
      fequ[2]=std::pow(var[2],2)/var[0]+pressure;
      fequ[3]=var[2]*var[3]/var[0];
      fequ[4]=(var[4]+pressure)*var[2]/var[0];    
    }
    
    void Euler_<3>::get_3d_fz(double* fequ, double* var){
      double pressure;
      pressure=get_thermal_pressure(var);        
      fequ[0]=var[3];
      fequ[1]=var[3]*var[1]/var[0];
      fequ[2]=var[3]*var[2]/var[0];
      fequ[3]=std::pow(var[3],2)/var[0]+pressure;
      fequ[4]=(var[4]+pressure)*var[3]/var[0];        
    }
    
  }
}
